Ab initio study of shock compressed oxygen 
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Quantum molecular dynamic simulations are introduced to study the shock compressed oxygen. 
The principal Hugoniot points derived from the equation of state agree well with the available 
experimental data. With the increase of pressure, molecular dissociation is observed. Electron 
spin polarization determines the electronic structure of the system under low pressure, while it is 
suppressed around 30 ~ 50 GPa. Particularly, nonmetal-metal transition is taken into account, 
which also occurs at about 30 ~ 50 GPa. In addition, the optical properties of shock compressed 
oxygen are also discussed. 

PACS numbers: 62.50.+P, 71.30.+h, 61.20.Ja 



I. INTRODUCTION 

The study of pressure-induced transformations of ma- 
terials in the shock compressed fluid state is of great in- 
terest in the last decade^. In particular, the homonu- 
clear diatomic molecular systems such as hydrogen 2-3 , 
oxygeni££i£, nitrogen^, and the halogens 10,11 , have 
been intensively studied with a variety of experimen- 
tal techniques, and the measurements of Hugoniots have 
reached megabar range. On the theoretical side, Quan- 
tum molecular dynamic (QMD) and path-integral Monte 
Carlo (PIMC) simulation s 12 ' 13 ! 14 have been introduced 
in the study of materials under extreme conditions. In 
spite of fact that these theoretical simulations in some 
cases already provided explanatory and predictive re- 
sults, to date, however, many fundamental questions of 
these above-mentioned systems in the fluid state are still 
under intense discussion due to the profound behavior 
they display. 

A knowledge of the properties of oxygen, in particu- 
lar, the equation of state (EOS) and derived properties, 
such as the Hugoniot curve and (D, u) diagram, are im- 
portant for shocks, detonations, biology, and fluids. Sin- 
gle or multiple shock-wave experiments have been per- 
formed on oxygen 4 ^, and principal Hugoniot are deter- 
mined from the EOS data. The partial molecular disso- 
ciation and a two-component conductive fluid are indi- 
cated by the results. Theoretically, a classical repulsive 
pair potential was introduced to investigate the EOS of 
shock-compressed oxygen^, and the dissociation effect 
of oxygen at high pressure was partially taken into ac- 
count by using a classical self-consistent fluid variational 
theory 1 -. Due to the intrinsic approximations of these 
classical methods, however, a full quantum-mechanical 
description of the change in the electronic structure of 
oxygen under high-pressure and high-temperature shock- 
wave compression still remains to be presented and un- 
derstood. The necessity for such quantum-mechanical 
treatment of shock compressed oxygen can, for exam- 
ple, be clearly seen by the following facts: (i) QMD 
method, where the electrons are fully quantum mechan- 



ical treated, has been proven successfully in calculat- 
ing the electronic structure, thermodynamical and opti- 
cal properties of warm dense matte r 17 ' 18 ; (ii) nonmetal- 
metal transition of oxygen was observed recently 6 ; (iii) 
low-pressure QMD calculations 1 ^ already indicated im- 
portance of electron spin polarization effect during oxy- 
gen dissociation. From these facts, therefore, elec- 
tronic structure should be seriously considered, the EOS 
data and Hugoniot points from density functional theory 
(DFT) are highly needed for shock compressed oxygen. 

In the present work, we apply QMD simulations to 
study the thermophysical properties of shock compressed 
oxygen. We determine the EOS data and dissociation 
fraction of oxygen along the principal Hugoniot. The pair 
correlation function is derived from the structural infor- 
mation during QMD simulations. The electronic struc- 
ture calculations within QMD provide the charge density 
in the simulation box and density of states (DOS) at ev- 
ery time step. The dissociation of oxygen and electronic 
structure are important for identifying nonmetal-metal 
transition. The dielectric function s(uj) and the conse- 
quent optical quantities are also extracted. 



II. QUANTUM MOLECULAR DYNAMICS 
SIMULATIONS 



In this study, we employ the VASP plane-wave pseu- 
dopotential code, which was developed at the Techni- 
cal University of Vienn a 19 ' 20 . The finite temperature 
density functional theory molecular dynamics (FTDFT- 
MD ) 21 ' 22 , where the electronic states are populated ac- 
cording to a Fermi-Dirac distribution at temperature T e , 
is used in the present work. The electron wave func- 
tions are calculated with the projector augmented wave 
(PAW) potentia l 23 ' 24 . The exchange correlation func- 
tional is calculated within generalized gradient approx- 
imation (GGA) using the parametrization of Perdew- 
Wang 91—. Electron spin polarization, which can be ef- 
fectively used to describe the character of O2 molecules, 
is also considered in our calculation. 
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In all the simulations, a total number of 96 oxygen 
atoms (48 2 molecules) are included in a supercell with 
periodic boundary condition. Trajectories are calculated 
at separate densities and temperatures. The selected 
densities range from p=2.5 to 3.6 g/cm 3 with temper- 
atures from 1000 to 16000 K. The plane wave cutoff en- 
ergy is set to be 600.0 eV. For molecular dynamic sim- 
ulations, only r point of the Brillouin zone is included, 
while 4x4x4 Monkhorst-Pack 2 ^ k points are used for 
calculating the electronic structure. Integration of equa- 
tions of motion proceeds 2000 steps with the time step of 
2 fs. Then, the system is equilibrated 200 steps and the 
final 300 steps are used to calculate physical quantities. 
During simulations, the ionic temperature T, is kept con- 
stant at every time step by using velocity scaling. The 
local thermodynamical equilibrium is reached by setting 
electron temperature T e equal to Tj. The accuracy of 
our calculations is examined by the bond length of O2 
molecule in its ground state, and the result is 1.24 A, 
which agrees well with the experimental. 



III. PRINCIPAL HUGONIOT AND PAIR 
CORRELATION FUNCTION 

A crucial measurement of the EOS data of oxygen un- 
der shock conditions is the Hugoniot^ which is the locus 
of points in (E, P, V)-space satisfying the condition: 



(Eq — Ei) 



\(Vq 



V 1 )(P +Pi) = 0, 



(1) 



where E is the internal energy, P is the pressure, V is 
the volume, and the subscripts and 1 denote the ini- 
tial and shocked state, respectively. This relation follows 
from conservation of matter, momentum, and energy for 
an isolated system compressed by a pusher at a constant 
velocity. In the canonical (NVT) ensemble in which both 
E and P are temperature dependent, the locus of states 
which satisfy Eq. (1) is the so-called principal Hugoniot, 
which describes the shock adiabat between the initial and 
final states. In our first-principles calculations, the in- 
ternal energy consists of the DFT total energy and the 
zero-point vibrational energy (\hv v n, per molecule). The 
pressure is evaluated via the forces provided by VASP. 
The particle velocity u p and shock velocity u s can be de- 
rived from the other two Rankinc-Hugoniot equations 28 
Vi = Vo[l — (up/u s )] and Pi — Pq = pqu s u p . In order to 
calculate the principal Hugoniot point for a given Vi, a 
series of simulations are executed for different tempera- 
tures T. {Eq-Ei) and | (V - V 1 ) (P + Pi ) are then fitted 
to polynomial expansion of T. The principal Hugoniot 
temperature Ti and pressure Pi are then determined by 
solving Eq. (1). 

As the starting point of the principal Hugoniot curve, 
it is important to determine the initial state (Eq, Pq, 
Vq). Generally, the initial-state internal energy should be 
calculated under the experimental condition, with p=1.2 
g/cm 3 and T=77 K. However, this calculation is highly 
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FIG. 1: (Color online). Left panel: Calculated principal 
Hugoniot curve, in comparison with the previous experimen- 
tal data; Right panel: u p -u s curve. 



TABLE I: Principal Hugoniot points derived from DFT-MD 
simulation at a series of density (p), pressure (P), and tem- 
perature (T). The corresponding particle velocity (u p ), shock 
velocity (u a ) and dissociation ratio (R) are included. 



P (g/cm 3 ) 


P (GPa) 


T(K) 


Up (km/s) 


Us (km/s) 


7? 


2.5 


19.4 


2413 


2.90 


5.57 


27% 


2.7 


32.8 


5258 


3.89 


7.01 


54% 


3.0 


49.8 


6489 


4.99 


8.31 


69% 


3.3 


71.9 


8823 


6.18 


9.70 


71% 



time demanding and is out of our computational capa- 
bilities. As an alternative selection, we determine the 
initial-state energy via single molecule energy. In the 
present work, we perform PAW potential calculation of 
an isolated (spin-polarized) O2 molecule at the equilib- 
rium bond length. With the inclusion of zero point en- 
ergy, the internal energy Eq is —4.84 eV/atom at T=77 
K. The initial pressure Pq can be neglected compared 
with the high pressures of shock compressed states. 

The principal Hugoniot curve and (u s , u p ) diagram are 
shown in Fig. 1. We find good agreement between our 
DFT-MD results and the previous experimental data^£. 
The present DFT-MD simulations can reasonably reflect 
the main tendency of oxygen thermodynamic properties 
along the principal Hugoniot curve. We now investigate 
explicit nature of the fluid for which the experiments pro- 
vide only indirect evidence. The probability of finding 
a atom at distance r from a reference atom could be 
given by pair correlation function g(r), which is obtained 
by averaging over all particles and simulation steps in 
equilibrium. The fraction of molecule can be derived by 
evaluating the coordination number, which could be ex- 
pressed as follows: 



K(r) 



N-l 
V 



Airr' 2 g(r')dr' . 



(2) 
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FIG. 3: (Color online). Isosurface of spin density for (a) p=2.5 
g/cm 3 , T=2400 K, (b) p=2.7 g/cm 3 , T=5200 K, (c) p=3.0 
g/cm 3 , T=6500 K, and (d) p=3.3 g/cm 3 , T=8800 K. 



FIG. 2: (Color online). Pair correlation functions (black line) 
and K(r) (red line) along the principal Hugoniot curve. 

The coordination number is a weighted integral over the 
pair correlation function g(r) of the ions. N is the to- 
tal number of ions and V is the volume of the supercell. 
The value K at the maximum of g(r) (around 1.24 A) is 
equal to the fraction of molecules in the supercell. The 
calculated Hugoniot points and related data, such as tem- 
peratures and dissociation fraction, are listed in Table I. 

Figure 2 shows the pair correlation functions and K(r) 
for four densities and temperatures. At the density ,0=2.5 
g/cm 3 and T=2400 K, the pair correlation function is 
featured by a sharp maximum peak around the equilib- 
rium bond length of 1.24 A, and K (1.24) = 73%. Under 
this condition, the pair correlation function exhibits the 
typical characteristics of a molecular fluid. As pressure 
increases, the first main peak begin to reduce in ampli- 
tude. This tendency means that the dissociation of O2 
molecules becomes remarkable. The molecule fraction 
is 46% at p=2.7 g/cm 3 and 31% at p=3.0 g/cm 3 . The 
crossover from molecular to atomic fluid exists around 
p=2.7 ~ 3.0 g/cm 3 , and the corresponding principal 
Hugoniot pressure range from 30 to 50 GPa. Under such 
conditions, O2 molecules dissociate gradually. 



IV. ELECTRONIC STRUCTURE 



gen at the four principal Hugoniot points. At low pres- 
sures (P = 20 ~ 30 GPa), the fluid mainly consists of 
disordered O2 molecules, and less than 50% fraction of 
molecules dissociate. In this case, the finite spin polar- 
ization characteristic of molecular oxygen is obvious, as 
shown in Fig. 3(a). As a result, the structural, electrical, 
and optical properties of the system is closely related with 
the electron spin polarization of oxygen molecules at low 
pressures. With the increase of pressures (P = 50 ~ 70 
GPa), electron spin polarization is suppressed, and the 
spin density in the supercell tends to disappear, as shown 
in Figs. 3(c) and 3(d). This change will consequently 
affect DOS of the system, which plays a key role in elu- 
cidating nonmetal-metal transition and optical proper- 
ties. As shown in Fig. 4, which plots spin-resolved DOS 
of shocked oxygen in different pressure regions, a large 
band gap exists in the DOS when the pressure is about 
20GPa (p=2.5 g/cm 3 , T=2400 K), and the system stays 
in the insulating state. With the increase of pressure, the 
electrons tend to be delocalized by the non-zero value of 
DOS at the Fermi energy. Finally, the localized molec- 
ular bonds as well as the valence-conduction band gap 
prominently disappear and the metallic like conductivity 
emerges at this stage. The nonmetal-metal transition, 
which is attributed to the dissociation of molecules and 
thermal broadening of the electronic states, occurs at the 
pressure range from 30 to 50 GPa. This result agrees well 
with the experiment 5 . 



Metallization of oxygen was first reported by Desgre- 
niers et al£^ in solid state, with the nonmetal-metal 
transition accomplished by structural change 30 . The 
nonmetal-metal transition of oxygen in fluid state was 
later reported by Bestea et al&. The ground state of oxy- 
gen is a spin-triplet state under low pressures. The molec- 
ular structure and magnetism properties mainly originate 
from electron spin polarization 31 . Therefore, the elec- 
tronic structure, which dominates the physical proper- 
ties, is examined along the principal Hugoniot curve. 

Figure 3 shows the isosurfaces of spin densities of oxy- 



V. OPTICAL PROPERTIES 

Having shown the EOS data and nonmetal-metal tran- 
sition along the principal Hugoniot curve, we turn now to 
study the optical properties of warm dense oxygen. The 
frequency dependent dielectric function e(w) has both 
real and imaginary parts: 

e(w) = e (1) M + e (2) H«. (3) 
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FIG. 4: Spin-resolved electron densities of states along the 
principal Hugoniot curve. The Fermi energy is set to be zero, 
(a) p=2.5g/cm 3 , T=2400 K; (b) p=2.7 g/cm 3 , T=5200 K; (c) 
p=3.0 g/cm 3 , T=6500 K; (d) p=3.3 g/cm 3 , T=8800 K. 



typical wavelength of 414 nm to investigate the change 
in reflectivity, as shown in the right panel of Fig. 5. 
One can see that the reflectivity increases from 0.01 to 
0.2 along the principal Hugoniot. The rapid increase in 
reflectivity occurs at the pressure range from 30 to 50 
GPa, with about 50% dissociation of molecules. 
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The real n(u>) and imaginary k(u>) parts of the refraction 
index are related to the dielectric function by 



e(w) = e (1) M + e (2) (iu)i = [n(u>) + k{u:)i\ 2 , 
or equivalenly 



l(w) = ^[| £ ( w )| + £ (i)H], 



k(u) 



|eH|-eW( W )]. 



(4) 



(5) 



(6) 



From the refraction index, the reflectivity r(uj) is given 
by 



r(w) = 



[1 - n(m)} 2 + kjiof 
[1 + n(uj)} 2 + jfc(w) 2 



(7) 



In Fig. 5 (left panel), we show the variation of the 
reflectivity r{oj) along the principal Hugoniot curve. At 
lower pressures, for example P=20^30 GPa, there are 
two peaks around the photon energies of 5.0 and 11.0 
eV, corresponding to the wavelength of 250 and 110 nm, 
respectively. With the increase in pressure, the peaks 
vanish. The reflectivity reaches a constant value of 0.172 
for the photon energy higher than 35.0eV. We select a 



FIG. 5: (Color online) Variation of the frequency-dependent 
optical reflectivity along the principal Hugoniots (left panel); 
Reflectivity at a fixed wavelength of 414 nm (right panel). 



VI. CONCLUSION 

In summary, we have applied QMD simulations in the 
study of shock-compressed oxygen. The principal Hugo- 
niot curve has been derived from the computational EOS 
data, which is consistent well with the experimental mea- 
surement. The electron spin polarization, determines 
physical properties of oxygen at lower pressures, has been 
shown to be greatly suppressed at higher pressures. This 
fact, as well as the nonmetal-metal transition and change 
in reflectivity found in this paper, are attributed to the 
dissociation of O2 molecules. We have determined that 
the molecular-atomic fluid transition exists at the pres- 
sure around 30^50 GPa, which also accord well with the 
recent experiment. 
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